Interpretation for Iq-tree concordance factors at http://www.robertlanfear.com/blog/files/archive-2018.html
First phylogenetic tree estimated with data from ortholog11-noBouvieri-queenconsambig.tar.gz, consensus alignments filtered through EAPhy, RAxML trees called as -s fasta -m GTRGAMMA -N 10 -p 12345 per locus and ASTRAL tree estimated (through script runRaxMLgeneTrees_ASTRALspeciesTree.R)
Showing tree comparison with published in Romiguier et al. 2017
tree <- read.newick("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allMessor_busco_iqtree/concord.cf.tree")
tree_root <- root(tree, outgroup = "SRR4292898", edgelabel = TRUE) ## by aciculatus as outgroup
nohyb <- dt %>%
filter(label %in% tree$tip.label)
grp <- list()
for (sp in unique(nohyb$species)) {
grp[[sp]] <- nohyb[nohyb$species %in% sp, "label"]
}
subgrp <- list()
for (sp in unique(nohyb$mitochondria)) {
subgrp[[sp]] <- MRCA(ggtree(tree_root), nohyb[nohyb$mitochondrial %in% sp, "label"])
}
subgrp <- subgrp[names(subgrp)[!names(subgrp) %in% ""]]
grouped <- groupOTU(tree_root, .node = grp)
p <- ggtree(grouped, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
nohyb[, c("label", "caste")] + geom_tree(aes(color = group)) + geom_tiplab(size = 2,
aes(color = group), hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
"#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, 28) +
geom_nodelab(size = 1.8, nudge_x = 0.18)
for (nd in 1:length(subgrp)) {
p <- p + geom_cladelabel(node = subgrp[[nd]], label = names(subgrp)[nd], offset = 2.2,
align = TRUE, fontsize = 2)
}
grid.arrange(p)iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allMessor_busco_iqtree/concord.cf.tree")
nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allMessor_busco_iqtree/geneTrees.treefile"))
pops3 <- left_join(data.frame(label = iqtree@phylo$tip.label), dt) %>%
mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
label2 = paste0(label, "_", lineage))
grp3 <- list()
for (sp in unique(pops3$species)) {
grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
.node = grp3)
treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops3[, c("label", "label2", "species", "caste")]
treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
size = 3, hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w"), utf8ToInt("m")),
breaks = c("queen", "worker", "male")) + scale_color_manual(values = c("#FFAD56",
"#C5EABD", "purple"), breaks = c("queen", "worker", "male")) + xlim(0, 25) +
geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))Using 2709 loci.
tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allMessor_busco_concatenatedRAxML/RAxML_bipartitionsBranchLabels.allMessor_busco')
tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)){
grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}
# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
# .node= grp)
grouped <- groupOTU(tree,.node= grp)
p <- ggtree(grouped, #branch.length='none',
ladderize=TRUE, right = TRUE)%<+%
dtTree[,c('label','caste','mitochondrial','label2','species')]
q1 <- p +
geom_tree(aes(color=group)) +
scale_color_manual(name = "species",
values = spCol[names(spCol) %in% names(grp)]) +
ggnewscale::new_scale_color() +
geom_tiplab(size=2, aes(label = label2, color = mitochondrial),
hjust = -0.09,
# align = TRUE,
linetype = "dotted", show.legend = F) +
scale_color_manual(values = linCol) +
guides(color=guide_legend(ncol=1)) +
ggnewscale::new_scale_color() +
geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) +
scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
breaks = c("queen", "male", "worker")) +
xlim(0, max(p$data$x)) +
geom_treescale() +
geom_label2(aes(label = bootstrap,
subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
theme(legend.position = "bottom")
q1tree <- read.raxml("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allNoHybrid_mitogenome_concatenatedRAxML/concatenated_raxml/RAxML_bipartitionsBranchLabels.mitogenome_75pConcatenated.tree")
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid, 1, 1), "_", label, "_", mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)) {
grp[[sp]] <- dtTree[dtTree$species %in% sp, "label"]
}
# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in%
# 'SRR4292897')), .node= grp)
grouped <- groupOTU(tree, .node = grp)
p <- ggtree(grouped, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
dtTree[, c("label", "caste", "mitochondrial", "label2", "species")]
q1 <- p + geom_tree(aes(color = group)) + scale_color_manual(name = "species", values = spCol[names(spCol) %in%
names(grp)]) + ggnewscale::new_scale_color() + geom_tiplab(size = 2, aes(label = label2,
color = mitochondrial), hjust = -0.09, aligned = TRUE, show.legend = F) + scale_color_manual(values = linCol) +
guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() + geom_tippoint(aes(x = x +
0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) + scale_shape_manual(values = c(utf8ToInt("q"),
utf8ToInt("m"), utf8ToInt("w")), breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
"#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, max(p$data$x) +
6) + geom_label2(aes(label = bootstrap, subset = as.numeric(sub("/.*", "", bootstrap)) >
50), size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) + theme(legend.position = "bottom")
q1tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allNoHybrid_busco_concatenatedRAxML/RAxML_bipartitionsBranchLabels.allNoHybrid_busco')
tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)){
grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}
# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
# .node= grp)
grouped <- groupOTU(tree,.node= grp)
p <- ggtree(grouped, #branch.length='none',
ladderize=TRUE, right = TRUE)%<+%
dtTree[,c('label','caste','mitochondrial','label2','species')]
q1 <- p +
geom_tree(aes(color=group)) +
scale_color_manual(name = "species",
values = spCol[names(spCol) %in% names(grp)]) +
ggnewscale::new_scale_color() +
geom_tiplab(size=2, aes(label = label2, color = mitochondrial),
hjust = -0.09,
# align = TRUE,
linetype = "dotted", show.legend = F) +
scale_color_manual(values = linCol) +
guides(color=guide_legend(ncol=1)) +
ggnewscale::new_scale_color() +
geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) +
scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
breaks = c("queen", "male", "worker")) +
xlim(0, max(p$data$x)) +
geom_treescale() +
geom_label2(aes(label = bootstrap,
subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
theme(legend.position = "bottom")
q1iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allNoHybrid_phasedH0_iqtree/iqtreeOutput_Filtered_concord/phasedH0_min53IndFiltered_concord.cf.tree")
pops3 <- dt %>%
mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
label = paste0(sra_access_number, ".0"), label2 = paste0(sra_access_number,
"_", lineage), label3 = sra_access_number) %>%
filter(label %in% iqtree@phylo$tip.label)
grp3 <- list()
for (sp in unique(pops3$species)) {
grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897.0")),
.node = grp3)
treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops3[, c("label", "label2", "species", "caste")]
treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
size = 3, hjust = -0.08, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
"#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, 21) +
geom_nodelab(size = 2, nudge_x = 0.42)tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/allNoMaleHybrid_phasedH0_concatenatedRAxML/RAxML_bipartitionsBranchLabels.allNoMaleHybrid_phasedH0.tree')
tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)){
grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}
# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
# .node= grp)
grouped <- groupOTU(tree,.node= grp)
p <- ggtree(grouped, #branch.length='none',
ladderize=TRUE, right = TRUE)%<+%
dtTree[,c('label','caste','mitochondrial','label2','species')]
q1 <- p +
geom_tree(aes(color=group)) +
scale_color_manual(name = "species",
values = spCol[names(spCol) %in% names(grp)]) +
ggnewscale::new_scale_color() +
geom_tiplab(size=2, aes(label = label2, color = mitochondrial),
hjust = -0.09,
# align = TRUE,
linetype = "dotted", show.legend = F) +
scale_color_manual(values = linCol) +
guides(color=guide_legend(ncol=1)) +
ggnewscale::new_scale_color() +
geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) +
scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
breaks = c("queen", "male", "worker")) +
xlim(0, max(p$data$x)) +
geom_treescale() +
geom_label2(aes(label = bootstrap,
subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
theme(legend.position = "bottom")
q1tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/AllNoMales_mitogenome_concatenatedRAxML/concatenated_raxml/RAxML_bipartitionsBranchLabels.mitogenome_75pConcatenated.tree')
tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)){
grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}
# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
# .node= grp)
grouped <- groupOTU(tree,.node= grp)
p <- ggtree(grouped, #branch.length='none',
ladderize=TRUE, right = TRUE)%<+%
dtTree[,c('label','caste','mitochondrial','label2','species')]
q1 <- p +
geom_tree(aes(color=group)) +
scale_color_manual(name = "species",
values = spCol[names(spCol) %in% names(grp)]) +
ggnewscale::new_scale_color() +
geom_tiplab(size=2, aes(label = label2, color = mitochondrial),
hjust = -0.09,
# align = TRUE,
linetype = "dotted", show.legend = F) +
scale_color_manual(values = linCol) +
guides(color=guide_legend(ncol=1)) +
ggnewscale::new_scale_color() +
geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) +
scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
breaks = c("queen", "male", "worker")) +
xlim(0, max(p$data$x)) +
geom_treescale() +
geom_label2(aes(label = bootstrap,
subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
theme(legend.position = "bottom")
q1tree <- read.newick("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/pureNoMales_mitogenome_iqtree/iqtreeOutput_concord/pureNoMales_concord.cf.tree")
tree_root <- root(tree, outgroup = "SRR4292898", edgelabel = TRUE) ## by aciculatus as outgroup
nohyb <- dt %>%
filter(sra_access_number %in% tree$tip.label)
grp <- list()
for (sp in unique(nohyb$species)) {
grp[[sp]] <- nohyb[nohyb$species %in% sp, "label"]
}
subgrp <- list()
for (sp in unique(nohyb$mitochondria)) {
subgrp[[sp]] <- MRCA(ggtree(tree_root), nohyb[nohyb$mitochondrial %in% sp, "label"])
}
subgrp <- subgrp[names(subgrp)[!names(subgrp) %in% ""]]
grouped <- groupOTU(tree_root, .node = grp)
p <- ggtree(grouped, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
nohyb[, c("label", "caste")] + geom_tree(aes(color = group)) + geom_tiplab(size = 2,
aes(color = group), hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
"#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, 20) +
geom_nodelab(size = 1.8, nudge_x = 0.18)
for (nd in 1:length(subgrp)) {
p <- p + geom_cladelabel(node = subgrp[[nd]], label = names(subgrp)[nd], offset = 2.2,
align = TRUE, fontsize = 2)
}
grid.arrange(p)tree <- read.iqtree('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/messorQueens_mitogenome_iqtree/concord.cf.tree')
nLoci <- nrow(read.delim('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_busco_iqtree/geneTrees.treefile'))
# tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)){
grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}
grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
.node= grp)
# tree2 <- tree
# tree2@phylo <- reorder(phangorn::midpoint(tree2@phylo)) ## to use a midpoint root
# grouped <- groupOTU(tree,.node= grp)
p <- ggtree(grouped, #branch.length='none',
ladderize=TRUE, right = TRUE)%<+%
dtTree[,c('label','caste','mitochondrial','label2','species')]
q1 <- p +
geom_tree(aes(color=group)) +
scale_color_manual(name = "species",
values = spCol[names(spCol) %in% names(grp)]) +
ggnewscale::new_scale_color() +
geom_tiplab(size=2, aes(label = label2, color = mitochondrial),
hjust = -0.09,
# align = TRUE,
linetype = "dotted", show.legend = F) +
scale_color_manual(values = linCol) +
guides(color=guide_legend(ncol=1)) +
ggnewscale::new_scale_color() +
xlim(0, max(p$data$x)+0.1) +
geom_treescale(x=0, y=0) +
geom_nodelab(aes(subset = as.numeric(sub("/.*", "", label)) > 50),
size = 2, hjust = 1.3, nudge_y = 0.15) +
theme(legend.position = "bottom")
# geom_label2(aes(label= label),subset = as.numeric(sub("/.*", "", label)) > 50)) + , size = 3, nudge_x = -0.5, label.padding = unit(0.1, "lines")
q1tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/messorQueens_mitogenome_concatenatedRAxML/RAxML_bipartitionsBranchLabels.messorQueens_mitogenome')
# tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)){
grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}
grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
.node= grp)
# grouped <- groupOTU(tree,.node= grp)
p <- ggtree(grouped, #branch.length='none',
ladderize=TRUE, right = TRUE)%<+%
dtTree[,c('label','caste','mitochondrial','label2','species')]
q1 <- p +
geom_tree(aes(color=group)) +
scale_color_manual(name = "species",
values = spCol[names(spCol) %in% names(grp)]) +
ggnewscale::new_scale_color() +
geom_tiplab(size=2, aes(label = label2, color = mitochondrial),
hjust = -0.09,
# align = TRUE,
linetype = "dotted", show.legend = F) +
scale_color_manual(values = linCol) +
guides(color=guide_legend(ncol=1)) +
xlim(0, max(p$data$x)+0.1) +
geom_treescale() +
geom_label2(aes(label = bootstrap,
subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
theme(legend.position = "bottom")
q1tree <- read.iqtree('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/messorQueens_busco_iqtree/concord.cf.tree')
nLoci <- nrow(read.delim('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_busco_iqtree/geneTrees.treefile'))
#tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)){
grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}
grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
.node= grp)
# grouped <- groupOTU(tree,.node= grp)
p <- ggtree(grouped, #branch.length='none',
ladderize=TRUE, right = TRUE)%<+%
dtTree[,c('label','caste','mitochondrial','label2','species')]
q1 <- p +
geom_tree(aes(color=group)) +
scale_color_manual(name = "species",
values = spCol[names(spCol) %in% names(grp)]) +
ggnewscale::new_scale_color() +
geom_tiplab(size=2, aes(label = label2, color = mitochondrial),
hjust = -0.09,
# align = TRUE,
linetype = "dotted", show.legend = F) +
scale_color_manual(values = linCol) +
guides(color=guide_legend(ncol=1)) +
xlim(0, max(p$data$x)) +
geom_treescale(x=0, y=0) +
geom_nodelab(aes(subset = as.numeric(sub("/.*", "", label)) > 50), size = 2,
hjust = 1.3, nudge_y = 0.15) +
theme(legend.position = "bottom")
# geom_label2(aes(label= label),subset = as.numeric(sub("/.*", "", label)) > 50)) + , size = 3, nudge_x = -0.5, label.padding = unit(0.1, "lines")
q1tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/messorQueens_busco_concatenatedRAxML/RAxML_bipartitionsBranchLabels.messorQueens_busco')
# tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)){
grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}
grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
.node= grp)
# grouped <- groupOTU(tree,.node= grp)
p <- ggtree(grouped, #branch.length='none',
ladderize=TRUE, right = TRUE)%<+%
dtTree[,c('label','caste','mitochondrial','label2','species')]
q1 <- p +
geom_tree(aes(color=group)) +
scale_color_manual(name = "species",
values = spCol[names(spCol) %in% names(grp)]) +
ggnewscale::new_scale_color() +
geom_tiplab(size=2, aes(label = label2, color = mitochondrial),
hjust = -0.09,
# align = TRUE,
linetype = "dotted", show.legend = F) +
scale_color_manual(values = linCol) +
guides(color=guide_legend(ncol=1)) +
xlim(0, max(p$data$x)) +
geom_treescale() +
geom_label2(aes(label = bootstrap,
subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
theme(legend.position = "bottom")
q1tree <- read.iqtree('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/messorQueens_phasedH0_iqtree/concord.cf.tree')
# tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)){
grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}
grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
.node= grp)
# grouped <- groupOTU(tree,.node= grp)
p <- ggtree(grouped, #branch.length='none',
ladderize=TRUE, right = TRUE)%<+%
dtTree[,c('label','caste','mitochondrial','label2','species')]
q1 <- p +
geom_tree(aes(color=group)) +
scale_color_manual(name = "species",
values = spCol[names(spCol) %in% names(grp)]) +
ggnewscale::new_scale_color() +
geom_tiplab(size=2, aes(label = label2, color = mitochondrial),
hjust = -0.09,
# align = TRUE,
linetype = "dotted", show.legend = F) +
scale_color_manual(values = linCol) +
guides(color=guide_legend(ncol=1)) +
xlim(0, max(p$data$x)+0.004) +
geom_treescale(x=0, y=0) +
geom_nodelab(aes(subset = as.numeric(sub("/.*", "", label)) > 50), size = 2, hjust = 1.3, nudge_y = 0.15) +
theme(legend.position = "bottom")
# geom_label2(aes(label= label),subset = as.numeric(sub("/.*", "", label)) > 50)) + , size = 3, nudge_x = -0.5, label.padding = unit(0.1, "lines")
q1tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/all/messorQueens_phasedH0_concatenatedRAxML/RAxML_bipartitionsBranchLabels.messorQueens_phasedH0')
# tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)){
grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}
grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
.node= grp)
# grouped <- groupOTU(tree,.node= grp)
p <- ggtree(grouped, #branch.length='none',
ladderize=TRUE, right = TRUE)%<+%
dtTree[,c('label','caste','mitochondrial','label2','species')]
q1 <- p +
geom_tree(aes(color=group)) +
scale_color_manual(name = "species",
values = spCol[names(spCol) %in% names(grp)]) +
ggnewscale::new_scale_color() +
geom_tiplab(size=2, aes(label = label2, color = mitochondrial),
hjust = -0.09,
# align = TRUE,
linetype = "dotted", show.legend = F) +
scale_color_manual(values = linCol) +
guides(color=guide_legend(ncol=1)) +
xlim(0, max(p$data$x)+0.01) +
geom_treescale(x=0, y=0) +
geom_nodelab(aes(subset = as.numeric(sub("/.*", "", label)) > 50), hjust = 1.3, nudge_y = 0.15, size = 2) +
theme(legend.position = "bottom")
# geom_label2(aes(label= label),subset = as.numeric(sub("/.*", "", label)) > 50)) + , size = 3, nudge_x = -0.5, label.padding = unit(0.1, "lines")
q1tree <- read.newick("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qsl_mitogenome_iqtree/iqtreeOutput_concord/qsl_concord.cf.tree")
tree_root <- midpoint(tree, node.labels = "support") ###midpoint root because for some reason forgot to include outgroup sample
nohyb <- dt %>%
filter(sra_access_number %in% tree$tip.label)
grp <- list()
for (sp in unique(nohyb$species)) {
grp[[sp]] <- nohyb[nohyb$species %in% sp, "label"]
}
subgrp <- list()
for (sp in unique(nohyb$mitochondria)) {
subgrp[[sp]] <- MRCA(ggtree(tree_root), nohyb[nohyb$mitochondrial %in% sp, "label"])
}
subgrp <- subgrp[names(subgrp)[!names(subgrp) %in% ""]]
grouped <- groupOTU(tree_root, .node = grp)
p <- ggtree(grouped, branch.length = "none", ladderize = TRUE, right = TRUE) + geom_tree(aes(color = group)) +
geom_tiplab(size = 3, aes(color = group), hjust = 0, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp)], name = "species") + xlim(0, 12) + geom_nodelab(size = 1.8, nudge_x = 0.5)
for (nd in 1:length(subgrp)) {
p <- p + geom_cladelabel(node = subgrp[[nd]], label = names(subgrp)[nd], offset = 1.7,
align = TRUE, fontsize = 3)
}
grid.arrange(p)path = "/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qsl_phased_AstralMrBayes/"
pops <- read.delim(paste0(path, "qsl_ID_pop.txt")) %>%
rename(label = ID) %>%
left_join(dt[, c("label", "species")])
tree <- read.tree(paste0(path, "/speciesTree/qslMrBayes_SpeciesTree.tree"))
grp <- list()
for (sp in unique(pops$species)) {
grp[[sp]] <- pops[pops$species %in% sp, "label"]
}
pops2 <- pops %>%
mutate(label2 = paste0(label, "_", lineage))
grouped <- groupOTU(root(tree, "SRR4292897"), .node = grp)
tree <- ggtree(grouped, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops2
tree$data$local.posterior <- ifelse(str_detect(as.numeric(tree$data$label), "[^a-zA-Z]"),
tree$data$label, "")
tree1 <- read.tree(paste0(path, "/speciesTree/qslMrBayes_SpeciesTree_t1.tree"))
grouped1 <- groupOTU(root(tree1, "SRR4292897"), .node = grp)
tree1 <- ggtree(grouped1, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops2
tree$data$quartet.support <- ifelse(str_detect(as.numeric(tree1$data$label), "[^a-zA-Z]"),
tree1$data$label, "")
tree8 <- read.astral(paste0(path, "/speciesTree/qslMrBayes_SpeciesTree_t8.tree"))
dt2 <- tree8@data %>%
select(q1:q3, node) %>%
drop_na()
pies <- nodepie(data.frame(dt2), cols = 1:3, color = c("red", "blue", "yellow"))
df <- tibble::tibble(node = as.numeric(dt2$node), pies = pies)
q <- tree %<+% df
library(ggpp)
q + geom_tiplab(aes(label = label2, color = group), size = 3, hjust = 0, aligned = TRUE) +
geom_plot(data = td_filter(!isTip), mapping = aes(x = x, y = y, label = pies),
vp.width = 0.1, vp.height = 0.07, hjust = 0.5, vjust = 0.5) + xlim(0, 12) +
scale_color_manual(values = spCol[names(spCol) %in% names(grp)]) + geom_text2(aes(label = as.numeric(local.posterior),
x = branch), size = 3, vjust = -0.2, color = "darkgreen") + geom_text2(aes(label = as.numeric(quartet.support),
x = branch), size = 3, vjust = 1.1, color = "darkred")iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qsl_phased_iqtree/phased_qsl.contree.tree")
iqtree <- treeio::drop.tip(iqtree, iqtree@phylo$tip.label[str_detect(iqtree@phylo$tip.label,
"\\.1")])
pops3 <- pops %>%
mutate(label2 = paste0(label, "_", lineage), label3 = label, label = paste0(label,
".0"))
grp3 <- list()
for (sp in unique(pops3$species)) {
grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897.0")),
.node = grp3)
treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops3
treeIQ + geom_tree(aes(color = group)) + geom_tiplab(aes(label = label2, color = group),
size = 3, hjust = 0, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + xlim(0, 13) + geom_nodelab(size = 3, nudge_x = 0.2)T1 <- tree + geom_tree(aes(color = group)) + geom_tiplab(aes(label = label2, color = group),
size = 2.5, hjust = 0, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + xlim(0, 15) + labs(title = "ASTRAL") + geom_nodelab(size = 3,
nudge_x = 0.12) + theme(legend.position = "none", plot.margin = unit(c(0, 0,
0, 0), "cm"))
d <- data.frame(old = pops3$label, new = pops3$label3)
groupedIQ@phylo = rename_taxa(groupedIQ@phylo, d, old, new)
pops32 <- pops3[, c("label3", "species", "label2")] %>%
rename(label = label3)
treeIQ2 <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops32
T2 <- treeIQ2 + geom_tree(aes(color = group)) + geom_tiplab(aes(label = label2, color = group),
size = 2.5, hjust = 1, aligned = TRUE) + scale_x_reverse(limits = c(13, 0)) +
scale_color_manual(values = spCol[names(spCol) %in% names(grp3)]) + labs(title = "IQ-tree") +
geom_nodelab(size = 3, nudge_x = 0.12) + theme(legend.position = "none", plot.margin = unit(c(0,
0, 0, 0), "cm"))
d1 = T1$data[T1$data$isTip, c(1:2, 4:10)]
d1$x[] = 1
d2 = T2$data[T2$data$isTip, c(1:4, 7:11)]
d2$x[] = 2
TTcon <- rbind(d1, d2)
L1 = ggplot(TTcon, aes(x = x, y = y, colour = group, group = label)) + geom_line() +
theme_void() + theme(legend.position = "none", plot.margin = unit(c(0, 0, 0,
0), "cm")) + scale_color_manual(values = spCol[names(spCol) %in% names(grp3)])
cowplot::plot_grid(T1, L1, T2, nrow = 1, align = "hv", rel_widths = c(1, 0.5, 1))qsl set of samples plus the three genome assemblies
iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslAssemblies_busco_iqtree/concord.cf.tree")
nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslAssemblies_busco_iqtree/geneTrees.treefile"))
pops3 <- left_join(data.frame(label = iqtree@phylo$tip.label), dt) %>%
mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
label2 = paste0(label, "_", lineage))
grp3 <- list()
for (sp in unique(pops3$species)) {
grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
.node = grp3)
treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops3[, c("label", "label2", "species", "caste")]
treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
size = 3, hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
"worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
"worker")) + xlim(0, 10) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))Using 2559 loci.
iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslAssemblies_phasedH0_iqtree/concord.cf.tree")
nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslAssemblies_busco_iqtree/geneTrees.treefile"))
pops3 <- left_join(data.frame(label = iqtree@phylo$tip.label), dt) %>%
mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
label2 = paste0(label, "_", lineage))
grp3 <- list()
for (sp in unique(pops3$species)) {
grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
.node = grp3)
treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops3[, c("label", "label2", "species", "caste")]
treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
size = 3, hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
"worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
"worker")) + xlim(0, 10) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))Using 2559 loci.
iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslAssemblies_phasedH01_iqtree/concord.cf.tree")
nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslAssemblies_phasedH01_iqtree/geneTrees.treefile"))
pops3 <- data.frame(ID = iqtree@phylo$tip.label) %>%
mutate(label = word(ID, sep = fixed("."), 1, 1)) %>%
left_join(dt) %>%
mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
label2 = paste0(ID, "_", lineage), label = ID)
grp3 <- list()
for (sp in unique(pops3$species)) {
grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
# groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in%
# 'SRR4292897')), .node= grp3)
groupedIQ <- groupOTU(iqtree, .node = grp3)
treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops3[, c("label", "label2", "species", "caste")]
treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
size = 3, hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
"worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
"worker")) + xlim(0, 18) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))Using 3543 loci.
qsl set of samples + a worker per lineage
iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslWorkerLineage_busco_iqtree/concord.cf.tree")
nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl/qslWorkerLineage_busco_iqtree/geneTrees.treefile"))
pops3 <- data.frame(ID = iqtree@phylo$tip.label) %>%
mutate(label = word(ID, sep = fixed("."), 1, 1)) %>%
left_join(dt) %>%
mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
label2 = paste0(ID, "_", lineage), label = ID)
grp3 <- list()
for (sp in unique(pops3$species)) {
grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
# groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in%
# 'SRR4292897')), .node= grp3)
groupedIQ <- groupOTU(iqtree, .node = grp3)
treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops3[, c("label", "label2", "species", "caste")]
treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
size = 3, hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
"worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
"worker")) + xlim(0, 18) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))Using 2640 loci.
Different samples chosen using only RNAseq dataset
Only one individual per species
beast_log_full <- parse_beast_tracelog_file("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2sub1sp_VCF_SNAPP/qsl2_vcfbedFiltered.log")
beast_log <- remove_burn_ins(beast_log_full, burn_in_fraction = 0.1)
sum_stats <- calc_summary_stats(beast_log, sample_interval = 1000)
kable(sum_stats, digits = 3)| mean | stderr_mean | stdev | variance | median | mode | geom_mean | hpd_interval_low | hpd_interval_high | act | ess | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| posterior | -10446.020 | 0.352 | 3.562 | 12.687 | -10445.608 | n/a | n/a | -10453.316 | -10439.850 | 17609.569 | 102.274 |
| likelihood | -10418.099 | 0.271 | 2.889 | 8.344 | -10417.808 | n/a | n/a | -10424.397 | -10413.449 | 15900.548 | 113.267 |
| prior | -27.921 | 0.191 | 2.032 | 4.127 | -27.543 | n/a | n/a | -32.237 | -24.675 | 15967.076 | 112.795 |
| lambda | 0.179 | 0.003 | 0.058 | 0.003 | 0.172 | n/a | 0.169961734922052 | 0.077 | 0.292 | 3782.016 | 476.201 |
| treeHeightLogger | 10.985 | 0.186 | 1.765 | 3.114 | 10.825 | n/a | 10.8479409107085 | 8.050 | 14.773 | 20065.133 | 89.758 |
| clockRate | 0.002 | 0.000 | 0.000 | 0.000 | 0.002 | n/a | 0.00176233759024145 | 0.001 | 0.002 | 22665.337 | 79.461 |
beast_log %>%
pivot_longer(cols = -Sample) %>%
ggplot(aes(value)) + geom_histogram() + facet_wrap(~name, scale = "free") + scale_x_continuous()beastMCC <- read.beast("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2sub1sp_VCF_SNAPP/qsl2_vcfbedFilteredMCC.tree")
p <- revts(ggtree(beastMCC, ladderize = TRUE))
p + geom_tiplab(aes(label = label), size = 4) + geom_range("height_0.95_HPD", color = "blue",
size = 1, alpha = 0.3) + geom_nodelab(aes(x = branch, label = round(posterior,
2), subset = posterior > 0.8), size = 4, size = 1.8, nudge_y = 0.2) + scale_x_continuous(labels = abs,
breaks = -(0:15)) + coord_cartesian(clip = "off") + theme_tree2(plot.margin = margin(6,
40, 6, 6))beast_log_full <- parse_beast_tracelog_file("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_VCF_SNAPPsp/qsl2_sp.log")
beast_log <- remove_burn_ins(beast_log_full, burn_in_fraction = 0.1)
sum_stats <- calc_summary_stats(beast_log, sample_interval = 1000)
kable(sum_stats, digits = 3)| mean | stderr_mean | stdev | variance | median | mode | geom_mean | hpd_interval_low | hpd_interval_high | act | ess | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| posterior | -7914.304 | 0.245 | 3.383 | 11.447 | -7914.002 | n/a | n/a | -7920.712 | -7907.756 | 18896.626 | 190.563 |
| likelihood | -7884.929 | 0.238 | 2.993 | 8.957 | -7884.626 | n/a | n/a | -7890.772 | -7879.766 | 22826.781 | 157.753 |
| prior | -29.375 | 0.066 | 1.516 | 2.297 | -29.067 | -28.4035981230664 | n/a | -32.592 | -27.080 | 6883.582 | 523.129 |
| lambda | 0.272 | 0.002 | 0.073 | 0.005 | 0.266 | 0.276314482283504 | 0.262524665428786 | 0.145 | 0.424 | 2064.477 | 1744.267 |
| treeHeightLogger | 7.108 | 0.046 | 0.700 | 0.490 | 7.068 | n/a | 7.07354539315322 | 5.789 | 8.527 | 15243.478 | 236.232 |
| clockRate | 0.004 | 0.000 | 0.000 | 0.000 | 0.004 | n/a | 0.00372456049255618 | 0.003 | 0.005 | 13701.272 | 262.822 |
beast_log %>%
pivot_longer(cols = -Sample) %>%
ggplot(aes(value)) + geom_histogram() + facet_wrap(~name, scale = "free") + scale_x_continuous()beastMCC <- read.beast("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_VCF_SNAPPsp/qsl2_sp_MCC.tree")
p <- revts(ggtree(beastMCC, ladderize = TRUE))
p + geom_tiplab(aes(label = label), size = 4) + geom_range("height_0.95_HPD", color = "blue",
size = 1, alpha = 0.3) + geom_nodelab(aes(x = branch, label = round(posterior,
2), subset = posterior > 0.8), size = 4, size = 1.8, nudge_y = 0.2) + scale_x_continuous(labels = abs,
breaks = -(0:15)) + coord_cartesian(clip = "off") + theme_tree2(plot.margin = margin(6,
40, 6, 6))beast_log_full <- parse_beast_tracelog_file("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_VCF_SNAPPsp-AciMono/run1/qsl2_sp_AciMono.log")
beast_log <- remove_burn_ins(beast_log_full, burn_in_fraction = 0.1)
sum_stats <- calc_summary_stats(beast_log, sample_interval = 1000)
kable(sum_stats, digits = 3)| mean | stderr_mean | stdev | variance | median | mode | geom_mean | hpd_interval_low | hpd_interval_high | act | ess | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| posterior | -16831.287 | 0.226 | 3.447 | 11.880 | -16830.815 | n/a | n/a | -16838.153 | -16825.185 | 3422.559 | 231.990 |
| likelihood | -16802.222 | 0.209 | 3.016 | 9.099 | -16801.902 | n/a | n/a | -16807.818 | -16796.998 | 3801.807 | 208.848 |
| prior | -29.064 | 0.102 | 1.723 | 2.970 | -28.671 | n/a | n/a | -32.521 | -26.495 | 2795.316 | 284.047 |
| lambda | 0.268 | 0.003 | 0.077 | 0.006 | 0.259 | n/a | 0.257946887799505 | 0.126 | 0.413 | 1460.433 | 543.674 |
| treeHeightLogger | 6.160 | 0.051 | 0.702 | 0.493 | 6.129 | n/a | 6.11966155006747 | 4.991 | 7.771 | 4110.355 | 193.171 |
| clockRate | 0.003 | 0.000 | 0.000 | 0.000 | 0.003 | n/a | 0.00332157547077767 | 0.003 | 0.004 | 4010.546 | 197.978 |
beast_log %>%
pivot_longer(cols = -Sample) %>%
ggplot(aes(value)) + geom_histogram() + facet_wrap(~name, scale = "free") + scale_x_continuous()beastMCC <- read.beast("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_VCF_SNAPPsp-AciMono/qsl2_sp_AciMonoCombinedMCC.tree")
p <- revts(ggtree(beastMCC, ladderize = TRUE))
p + geom_tiplab(aes(label = label), size = 4) + geom_range("height_0.95_HPD", color = "blue",
size = 1, alpha = 0.3) + geom_nodelab(aes(x = branch, label = round(posterior,
2), subset = posterior > 0.8), size = 4, size = 1.8, nudge_y = 0.2) + scale_x_continuous(labels = abs,
breaks = -(0:15)) + coord_cartesian(clip = "off") + theme_tree2(plot.margin = margin(6,
40, 6, 6))beast_log_full <- parse_beast_tracelog_file("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2-noAci_VCF_SNAPPsp/qsl2_sp_noAci.log")
beast_log <- remove_burn_ins(beast_log_full, burn_in_fraction = 0.1)
sum_stats <- calc_summary_stats(beast_log, sample_interval = 1000)
kable(sum_stats, digits = 3)| mean | stderr_mean | stdev | variance | median | mode | geom_mean | hpd_interval_low | hpd_interval_high | act | ess | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| posterior | -8099.867 | 0.139 | 3.117 | 9.718 | -8099.446 | n/a | n/a | -8106.267 | -8094.486 | 8969.918 | 501.788 |
| likelihood | -8071.497 | 0.140 | 2.860 | 8.179 | -8071.180 | n/a | n/a | -8077.330 | -8066.469 | 10722.871 | 419.757 |
| prior | -28.370 | 0.074 | 1.774 | 3.149 | -28.123 | n/a | n/a | -32.168 | -25.422 | 7792.797 | 577.585 |
| lambda | 0.277 | 0.002 | 0.080 | 0.006 | 0.269 | 0.279568483126404 | 0.265861374445471 | 0.139 | 0.441 | 1676.720 | 2684.408 |
| treeHeightLogger | 6.980 | 0.034 | 0.684 | 0.467 | 6.931 | n/a | 6.94737794002518 | 5.626 | 8.284 | 11231.668 | 400.742 |
| clockRate | 0.003 | 0.000 | 0.000 | 0.000 | 0.003 | n/a | 0.00300838695349966 | 0.002 | 0.004 | 7267.469 | 619.335 |
beast_log %>%
pivot_longer(cols = -Sample) %>%
ggplot(aes(value)) + geom_histogram() + facet_wrap(~name, scale = "free") + scale_x_continuous()beastMCC <- read.beast("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2-noAci_VCF_SNAPPsp/qsl2_sp_noAci_MCC.tree")
p <- revts(ggtree(beastMCC, ladderize = TRUE))
p + geom_tiplab(aes(label = label), size = 4) + geom_range("height_0.95_HPD", color = "blue",
size = 1, alpha = 0.3) + geom_nodelab(aes(x = branch, label = round(posterior,
2), subset = posterior > 0.8), size = 4, size = 1.8, nudge_y = 0.2) + scale_x_continuous(labels = abs,
breaks = -(0:15)) + coord_cartesian(clip = "off") + theme_tree2(plot.margin = margin(6,
40, 6, 6))beast_log_full <- parse_beast_tracelog_file("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2-noAciAre_VCF_spSNAPP/run1/qsl2_sp_noAciAre.log")
beast_log <- remove_burn_ins(beast_log_full, burn_in_fraction = 0.1)
sum_stats <- calc_summary_stats(beast_log, sample_interval = 1000)
kable(sum_stats, digits = 3)| mean | stderr_mean | stdev | variance | median | mode | geom_mean | hpd_interval_low | hpd_interval_high | act | ess | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| posterior | -16039.148 | 0.154 | 3.218 | 10.354 | -16038.733 | n/a | n/a | -16045.730 | -16033.573 | 10264.242 | 438.513 |
| likelihood | -16013.579 | 0.135 | 2.805 | 7.866 | -16013.271 | n/a | n/a | -16019.087 | -16008.492 | 10501.405 | 428.609 |
| prior | -25.569 | 0.074 | 1.562 | 2.440 | -25.216 | n/a | n/a | -28.756 | -23.214 | 10214.810 | 440.635 |
| lambda | 0.322 | 0.002 | 0.090 | 0.008 | 0.313 | n/a | 0.310029883926007 | 0.165 | 0.506 | 1668.415 | 2697.771 |
| treeHeightLogger | 6.854 | 0.032 | 0.618 | 0.382 | 6.828 | 7.58032501811522 | 6.82618238315202 | 5.721 | 8.119 | 11876.655 | 378.979 |
| clockRate | 0.005 | 0.000 | 0.000 | 0.000 | 0.005 | n/a | 0.00475405329388956 | 0.004 | 0.006 | 12584.100 | 357.674 |
beast_log %>%
pivot_longer(cols = -Sample) %>%
ggplot(aes(value)) + geom_histogram() + facet_wrap(~name, scale = "free") + scale_x_continuous()beastMCC <- read.beast("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2-noAciAre_VCF_spSNAPP/qsl2_sp_noAciAreCombinedMCC.tree")
p <- revts(ggtree(beastMCC, ladderize = TRUE))
p + geom_tiplab(aes(label = label), size = 4) + geom_range("height_0.95_HPD", color = "blue",
size = 1, alpha = 0.3) + geom_nodelab(aes(x = branch, label = round(posterior,
2), subset = posterior > 0.8), size = 4, size = 1.8, nudge_y = 0.2) + scale_x_continuous(labels = abs,
breaks = -(0:15)) + coord_cartesian(clip = "off") + theme_tree2(plot.margin = margin(6,
40, 6, 6))beast_log_full <- parse_beast_tracelog_file("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_VCF_SNAPPsample/qsl2_sample.log")
beast_log <- remove_burn_ins(beast_log_full, burn_in_fraction = 0.1)
sum_stats <- calc_summary_stats(beast_log, sample_interval = 1000)
kable(sum_stats, digits = 3)| mean | stderr_mean | stdev | variance | median | mode | geom_mean | hpd_interval_low | hpd_interval_high | act | ess | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| posterior | -14875.395 | 0.480 | 5.584 | 31.181 | -14874.934 | n/a | n/a | -14886.017 | -14864.935 | 18999.297 | 135.110 |
| likelihood | -14825.602 | 0.469 | 4.799 | 23.026 | -14825.243 | n/a | n/a | -14834.713 | -14817.108 | 24508.675 | 104.738 |
| prior | -49.794 | 0.108 | 2.843 | 8.084 | -49.388 | n/a | n/a | -55.509 | -45.099 | 3722.448 | 689.600 |
| lambda | 0.518 | 0.002 | 0.101 | 0.010 | 0.510 | n/a | 0.508478389893188 | 0.336 | 0.728 | 1355.567 | 1893.672 |
| treeHeightLogger | 8.200 | 0.026 | 0.766 | 0.587 | 8.163 | n/a | 8.16470357121235 | 6.667 | 9.600 | 3031.414 | 846.800 |
| clockRate | 0.003 | 0.000 | 0.000 | 0.000 | 0.003 | n/a | 0.0034195990990794 | 0.003 | 0.004 | 2850.903 | 900.416 |
beast_log %>%
pivot_longer(cols = -Sample) %>%
ggplot(aes(value)) + geom_histogram() + facet_wrap(~name, scale = "free") + scale_x_continuous()beastMCC <- read.beast("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_VCF_SNAPPsample/qsl2_sampleMCC.tree")
p <- revts(ggtree(beastMCC, ladderize = TRUE))
p + geom_tiplab(aes(label = label), size = 4) + geom_range("height_0.95_HPD", color = "blue",
size = 1, alpha = 0.3) + geom_nodelab(aes(x = branch, label = round(posterior,
2), subset = posterior > 0.8), size = 4, nudge_y = 0.2) + scale_x_continuous(labels = abs,
breaks = -(0:15)) + coord_cartesian(clip = "off") + theme_tree2(plot.margin = margin(6,
120, 6, 6))tree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_mitogenome_iqtree/concord.cf.tree")
nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_mitogenome_iqtree/geneTrees.treefile"))
pops3 <- dt %>%
mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
label2 = paste0(sra_access_number, "_", lineage), label = sra_access_number) %>%
filter(label %in% iqtree@phylo$tip.label)
grp3 <- list()
for (sp in unique(pops3$species)) {
grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
.node = grp3)
treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops3[, c("label", "label2", "species", "caste")]
treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
size = 3, hjust = -0.08, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.11, shape = caste, color = caste), size = 3, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
"worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
"worker")) + xlim(0, 15) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_busco_concatenatedRAxML/RAxML_bipartitionsBranchLabels.qsl2_busco_concatenated.tree')
tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)){
grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}
# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
# .node= grp)
grouped <- groupOTU(tree,.node= grp)
p <- ggtree(grouped, #branch.length='none',
ladderize=TRUE, right = TRUE)%<+%
dtTree[,c('label','caste','mitochondrial','label2','species')]
q1 <- p +
geom_tree(aes(color=group)) +
scale_color_manual(name = "species",
values = spCol[names(spCol) %in% names(grp)]) +
ggnewscale::new_scale_color() +
geom_tiplab(size=2, aes(label = label2, color = mitochondrial),
hjust = -0.09,
# align = TRUE,
linetype = "dotted", show.legend = F) +
scale_color_manual(values = linCol) +
guides(color=guide_legend(ncol=1)) +
ggnewscale::new_scale_color() +
geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) +
scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
breaks = c("queen", "male", "worker")) +
xlim(0, max(p$data$x)) +
geom_treescale() +
geom_label2(aes(label = bootstrap,
subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
theme(legend.position = "bottom")
q1tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_mitogenome_concatenatedRAxML/RAxML_bipartitionsBranchLabels.qsl2_mitogenome')
tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)){
grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}
# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
# .node= grp)
grouped <- groupOTU(tree,.node= grp)
p <- ggtree(grouped, #branch.length='none',
ladderize=TRUE, right = TRUE)%<+%
dtTree[,c('label','caste','mitochondrial','label2','species')]
q1 <- p +
geom_tree(aes(color=group)) +
scale_color_manual(name = "species",
values = spCol[names(spCol) %in% names(grp)]) +
ggnewscale::new_scale_color() +
geom_tiplab(size=2, aes(label = label2, color = mitochondrial),
hjust = -0.09,
# align = TRUE,
linetype = "dotted", show.legend = F) +
scale_color_manual(values = linCol) +
guides(color=guide_legend(ncol=1)) +
ggnewscale::new_scale_color() +
geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) +
scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
breaks = c("queen", "male", "worker")) +
xlim(0, max(p$data$x)) +
geom_treescale() +
geom_label2(aes(label = bootstrap,
subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
theme(legend.position = "bottom")
q1tree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_busco_iqtree/concord.cf.tree")
nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_busco_iqtree/geneTrees.treefile"))
pops3 <- dt %>%
mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
label2 = paste0(sra_access_number, "_", lineage), label = sra_access_number) %>%
filter(label %in% iqtree@phylo$tip.label)
grp3 <- list()
for (sp in unique(pops3$species)) {
grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
.node = grp3)
treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops3[, c("label", "label2", "species", "caste")]
treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
size = 3, hjust = -0.08, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.11, shape = caste, color = caste), size = 3, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
"worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
"worker")) + xlim(0, 15) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2_phasedH0_concatenatedRAxML/RAxML_bipartitionsBranchLabels.qsl2RNAseq_phasedH0')
tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)){
grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}
# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
# .node= grp)
grouped <- groupOTU(tree,.node= grp)
p <- ggtree(grouped, #branch.length='none',
ladderize=TRUE, right = TRUE)%<+%
dtTree[,c('label','caste','mitochondrial','label2','species')]
q1 <- p +
geom_tree(aes(color=group)) +
scale_color_manual(name = "species",
values = spCol[names(spCol) %in% names(grp)]) +
ggnewscale::new_scale_color() +
geom_tiplab(size=2, aes(label = label2, color = mitochondrial),
hjust = -0.09,
# align = TRUE,
linetype = "dotted", show.legend = F) +
scale_color_manual(values = linCol) +
guides(color=guide_legend(ncol=1)) +
ggnewscale::new_scale_color() +
geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) +
scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
breaks = c("queen", "male", "worker")) +
xlim(0, max(p$data$x)) +
geom_treescale() +
geom_label2(aes(label = bootstrap,
subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
theme(legend.position = "bottom")
q1With two samples as outgroup
tree <- read.raxml('/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/qsl2/qsl2-2out_phasedH0_concatenatedRAxML/RAxML_bipartitionsBranchLabels.qsl2-2out_phasedH0')
tree <- treeio::drop.tip(tree, c("SRR4292897","SRR4292898"))
dtTree <- dt %>%
filter(label %in% tree@phylo$tip.label) %>%
mutate(label2 = paste0(substring(hybrid,1,1),"_", label, "_",mitochondrial))
grp <- list()
for (sp in unique(dtTree$species)){
grp[[sp]] <- dtTree[dtTree$species %in% sp,'label']
}
# grouped <- groupOTU(root(tree, which(tree@phylo$tip.label %in% "SRR4292897")),
# .node= grp)
grouped <- groupOTU(tree,.node= grp)
p <- ggtree(grouped, #branch.length='none',
ladderize=TRUE, right = TRUE)%<+%
dtTree[,c('label','caste','mitochondrial','label2','species')]
q1 <- p +
geom_tree(aes(color=group)) +
scale_color_manual(name = "species",
values = spCol[names(spCol) %in% names(grp)]) +
ggnewscale::new_scale_color() +
geom_tiplab(size=2, aes(label = label2, color = mitochondrial),
hjust = -0.09,
# align = TRUE,
linetype = "dotted", show.legend = F) +
scale_color_manual(values = linCol) +
guides(color=guide_legend(ncol=1)) +
ggnewscale::new_scale_color() +
geom_tippoint(aes(x=x+0.12, shape = caste, color = caste),
size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) +
scale_color_manual(values = c( "#FFAD56" ,"#59A3F7", "#C5EABD"),
breaks = c("queen", "male", "worker")) +
xlim(0, max(p$data$x)) +
geom_treescale() +
geom_label2(aes(label = bootstrap,
subset = as.numeric(sub("/.*", "", bootstrap)) > 50),
size = 3, nudge_x = 0, label.padding = unit(0.1, "lines")) +
theme(legend.position = "bottom")
q1iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MbarDecCapOut_busco_iqtree/concord.cf.tree")
nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MbarDecCapOut_busco_iqtree/geneTrees.treefile"))
pops3 <- dt %>%
mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
label2 = paste0(sra_access_number, "_", lineage), label = sra_access_number) %>%
filter(label %in% iqtree@phylo$tip.label)
grp3 <- list()
for (sp in unique(pops3$species)) {
grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
.node = grp3)
treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops3[, c("label", "label2", "species", "caste")]
treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
size = 3, hjust = -0.08, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.11, shape = caste, color = caste), size = 3, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
"worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
"worker")) + xlim(0, 15) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))Using 2328 loci.
iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MbarQMWorkerLineage_busco_iqtree/concord.cf.tree")
nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MbarQMWorkerLineage_busco_iqtree/geneTrees.treefile"))
pops3 <- dt %>%
mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
label2 = paste0(sra_access_number, "_", lineage), label = sra_access_number) %>%
filter(label %in% iqtree@phylo$tip.label)
grp3 <- list()
for (sp in unique(pops3$species)) {
grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
.node = grp3)
treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops3[, c("label", "label2", "species", "caste")]
treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
size = 3, hjust = -0.08, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.11, shape = caste, color = caste), size = 3, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
"#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, 11) +
geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))Using 2263 loci.
iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MebeQ123Out_busco_iqtree/concord.cf.tree")
nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MebeQ123Out_busco_iqtree/geneTrees.treefile"))
pops3 <- dt %>%
mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
label2 = paste0(sra_access_number, "_", lineage), label = sra_access_number) %>%
filter(label %in% iqtree@phylo$tip.label)
grp3 <- list()
for (sp in unique(pops3$species)) {
grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
.node = grp3)
treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops3[, c("label", "label2", "species", "caste")]
treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
size = 3, hjust = -0.08, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.035, shape = caste, color = caste), size = 3, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
"worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
"worker")) + xlim(0, 5) + geom_nodelab(size = 3, nudge_x = 0.3) + guides(color = guide_legend(ncol = 1))Using 1654 loci.
iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MebeQ123Out_busco_iqtree/concord.cf.tree")
nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MebeQ123Out_busco_iqtree/geneTrees.treefile"))
pops3 <- dt %>%
mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
label2 = paste0(sra_access_number, "_", lineage), label = sra_access_number) %>%
filter(label %in% iqtree@phylo$tip.label)
grp3 <- list()
for (sp in unique(pops3$species)) {
grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
.node = grp3)
treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops3[, c("label", "label2", "species", "caste")]
treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
size = 3, hjust = -0.09, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("w")), breaks = c("queen",
"worker")) + scale_color_manual(values = c("#FFAD56", "#C5EABD"), breaks = c("queen",
"worker")) + xlim(0, 5) + geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))Using 1654 loci.
iqtree <- read.iqtree("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MebeOut_amb_iqtree/concord.cf.tree")
nLoci <- nrow(read.delim("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/MebeOut_amb_iqtree/geneTrees.treefile"))
pops3 <- dt %>%
mutate(lineage = case_when(mitochondrial %in% "" ~ species, TRUE ~ mitochondrial),
label2 = paste0(sra_access_number, "_", lineage), label = sra_access_number) %>%
filter(label %in% iqtree@phylo$tip.label)
grp3 <- list()
for (sp in unique(pops3$species)) {
grp3[[sp]] <- pops3[pops3$species %in% sp, "label"]
}
groupedIQ <- groupOTU(root(iqtree, which(iqtree@phylo$tip.label %in% "SRR4292897")),
.node = grp3)
treeIQ <- ggtree(groupedIQ, branch.length = "none", ladderize = TRUE, right = TRUE) %<+%
pops3[, c("label", "label2", "species", "caste")]
treeIQ + geom_tree(aes(color = species)) + geom_tiplab(aes(label = label2, color = species),
size = 3, hjust = -0.08, aligned = TRUE) + scale_color_manual(values = spCol[names(spCol) %in%
names(grp3)]) + guides(color = guide_legend(ncol = 1)) + ggnewscale::new_scale_color() +
geom_tippoint(aes(x = x + 0.11, shape = caste, color = caste), size = 3, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
"#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, 13) +
geom_nodelab(size = 2, nudge_x = 0.42) + guides(color = guide_legend(ncol = 1))Using 1601 loci.
tree <- read.newick("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/allIbericus_mitogenome_iqtree/concord.cf.tree")
nohyb <- dt %>%
filter(sra_access_number %in% tree$tip.label)
p <- ggtree(tree, branch.length = "none", ladderize = TRUE, right = TRUE) %<+% nohyb[,
c("label", "caste")] + geom_tree() + geom_tiplab(size = 2, hjust = -0.09, aligned = TRUE) +
geom_tippoint(aes(x = x + 0.12, shape = caste, color = caste), size = 2, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
"#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, 25) +
geom_nodelab(size = 1.8, nudge_x = 0.4)
grid.arrange(p)tree <- read.newick("/Users/acas/Dropbox/Post-docs/Roux_Lab/messor_analyses/phylo_analyses/analyses/perClade/ibericusQueens_phasedH0_iqtree/concord.cf.tree")
nohyb <- dt %>%
filter(sra_access_number %in% tree$tip.label)
p <- ggtree(tree, branch.length = "none", ladderize = TRUE, right = TRUE) %<+% nohyb[,
c("label", "caste")] + geom_tree() + geom_tiplab(size = 3, hjust = -0.09, aligned = TRUE) +
geom_tippoint(aes(x = x + 0.05, shape = caste, color = caste), size = 3, na.rm = TRUE) +
scale_shape_manual(values = c(utf8ToInt("q"), utf8ToInt("m"), utf8ToInt("w")),
breaks = c("queen", "male", "worker")) + scale_color_manual(values = c("#FFAD56",
"#59A3F7", "#C5EABD"), breaks = c("queen", "male", "worker")) + xlim(0, 10) +
geom_nodelab(size = 1.8, nudge_x = 0.4)
grid.arrange(p)